clc; clear; close all;

% 参数设置
N_iter = 500;         % 总迭代步数
N_discard = 300;       % 丢弃初始过渡序列
m = 2;                 % 样本熵嵌入维度

r_list = linspace(3, 4, 50);   % 参数 r 范围
b_list = linspace(0, 30, 50);  % 参数 b 范围

SE_mat = zeros(length(r_list), length(b_list)); % 存储样本熵结果

% 初始复数状态
% x0 = 0.463442265 + 0i;
% y0 = 0.04532285  + 0i;
% z0 = 0.002136285 + 0i;
x0 = 0.3;
y0 = 0.5;
z0 = 0.6;

% CTBCS固定参数
r1 = -10;
b1 = -3;
a = 1;
b2 = 2*pi;
c = 11.5;

parfor ir = 1:length(r_list)
    r = r_list(ir);
    SE_tmp = zeros(1, length(b_list));
    for ib = 1:length(b_list)
        b = b_list(ib);

        x = x0; y = y0; z = z0;
        X = zeros(N_iter,1);

        for k = 1:N_iter
            state = CTBCS(x, y, z, r, b, r1, b1, a, b2, c);
            x = state(1);
            y = state(2);
            z = state(3);

            X(k) = real(x); % 取实部计算样本熵
        end

        data = X(N_discard+1:end);
        r_tol = 0.2 * std(data);
        SE_tmp(ib) = SE(data, m, r_tol); % 调用你们封装的样本熵函数
    end
    SE_mat(ir, :) = SE_tmp;
end

% 归一化到0-1
SE_mat = (SE_mat - min(SE_mat(:))) / (max(SE_mat(:)) - min(SE_mat(:)));

% 网格坐标
[B_grid, R_grid] = meshgrid(b_list, r_list);

% 绘制三维曲面图
figure;
surf(B_grid, R_grid, SE_mat', 'EdgeColor', 'none');
set(gca, 'YDir', 'reverse'); 
xlabel('r');
ylabel('b');
zlabel('SE');
title('CTBCS 系统的 SE 复杂度');
colorbar;
view(45, 30);
grid on;

